load objmat;
trisurf(S.face, S.x , S.y ,S.z); hold on 
axis equal
    exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\tumor.obj', [S.x S.y S.z], [S.face]);
if(strcmp(S.method,'surrogate'))
    load Results.mat
    v=Data.xbest;        % surrogate model
    fbest=Data.fbest;
elseif(strcmp(S.method,'ES'))
    v= min_x{1,gen}(:,1); % ES
    fbest = min_f{1,gen}(1,1);
elseif(strcmp(S.method,'sim'))
    v=x0;                % sim annealing
    fbest=f0;
elseif(strcmp(S.method,'MCS'))
    v=pg;
    fbest=F{S.runs,1}(1,1);
elseif(strcmp(S.method,'Pattern'))
    v=x;
    fbest=fval;
elseif(strcmp(S.method,'GA'))
    v=x;
    fbest=fval;
elseif(strcmp(S.method,'simann'))
    v=x;
    fbest=fval;
end
xr=S.xr;
%%
if(S.bioheat==0)
    [x1,y1,z1]=ellipsoid(v(1),v(2),v(3),xr(1),xr(2),xr(3));
    xEllipsoid1 = reshape(x1,1,[])';
    yEllipsoid1 = reshape(y1,1,[])';
    zEllipsoid1 = reshape(z1,1,[])';
    teta1=deg2rad(v(7));
    phi1=deg2rad(v(8));
    xrotEllipsoid1 = [xEllipsoid1-v(1) yEllipsoid1-v(2) zEllipsoid1-v(3)]*[1 0 0;0 cos(teta1) -sin(teta1);0 sin(teta1) cos(teta1)];
    yrotEllipsoid1 = xrotEllipsoid1*[cos(phi1) 0 sin(phi1);0 1 0;-sin(phi1) 0 cos(phi1)];
    face1=delaunay(yrotEllipsoid1);
    title(fbest, 'fontsize', 10);
    trisurf(face1,yrotEllipsoid1(:,1)+v(1),yrotEllipsoid1(:,2)+v(2),yrotEllipsoid1(:,3)+v(3),'FaceColor','flat'); hold on
    axis equal
%     exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\iceball_1.obj', [yrotEllipsoid1(:,1)+v(1),yrotEllipsoid1(:,2)+v(2),yrotEllipsoid1(:,3)+v(3)],[])
    %%
    if(S.needle_num >= 2)
    [x2,y2,z2]=ellipsoid(v(4),v(5),v(6),xr(1),xr(2),xr(3));
    teta2=deg2rad(v(9));
    phi2=deg2rad(v(10));
    xEllipsoid2 = reshape(x2,1,[])';
    yEllipsoid2 = reshape(y2,1,[])';
    zEllipsoid2 = reshape(z2,1,[])';
    xrotEllipsoid2 = [xEllipsoid2-v(4) yEllipsoid2-v(5) zEllipsoid2-v(6)]*[1 0 0;0 cos(teta2) -sin(teta2);0 sin(teta2) cos(teta2)];
    yrotEllipsoid2 = xrotEllipsoid2*[cos(phi2) 0 sin(phi2);0 1 0;-sin(phi2) 0 cos(phi2)];
    face2=delaunay(yrotEllipsoid2);
    trisurf(face2,yrotEllipsoid2(:,1)+v(4),yrotEllipsoid2(:,2)+v(5),yrotEllipsoid2(:,3)+v(6),'FaceColor','flat'); 
    axis equal
%         exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\iceball_2.obj', [yrotEllipsoid2(:,1)+v(4),yrotEllipsoid2(:,2)+v(5),yrotEllipsoid2(:,3)+v(6)], [])
    end
    %%
    if(S.needle_num >= 3)
        [x3,y3,z3]=ellipsoid(v(11),v(12),v(13),xr(1),xr(2),xr(3));
        teta3=deg2rad(v(14));
        phi3=deg2rad(v(15));
        xEllipsoid3 = reshape(x3,1,[])';
        yEllipsoid3 = reshape(y3,1,[])';
        zEllipsoid3 = reshape(z3,1,[])';
        xrotEllipsoid3 = [xEllipsoid3-v(11) yEllipsoid3-v(12) zEllipsoid3-v(13)]*[1 0 0;0 cos(teta3) -sin(teta3);0 sin(teta3) cos(teta3)];
        yrotEllipsoid3 = xrotEllipsoid3*[cos(phi3) 0 sin(phi3);0 1 0;-sin(phi3) 0 cos(phi3)];
        face3=delaunay(yrotEllipsoid3);
        trisurf(face3,yrotEllipsoid3(:,1)+v(11),yrotEllipsoid3(:,2)+v(12),yrotEllipsoid3(:,3)+v(13),'FaceColor','flat'); 
        axis equal
%         exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\iceball_3.obj', [yrotEllipsoid3(:,1)+v(11),yrotEllipsoid3(:,2)+v(12),yrotEllipsoid3(:,3)+v(13)],[])
    end
    hold off
end
%%
if(S.bioheat==1)
    if(isempty(S.vessels))
        b=mexcpp(v,S.time,S.needle_num,S.xr);
    else
        b=mexcpp(v,S.time,S.needle_num,S.xr,S.row_vessel-1,S.col_vessel-1,S.I3_vessel-1);
    end
    BW=b<S.temp;
    in=0;
    for i=1:size(S.row,1)
        B=sum(BW(S.row(i,1)+S.boundx,S.col(i,1)+S.boundy,S.I3(i,1)+S.boundz));
        in=in+B;
    end
    A=sum(sum(sum(BW)));
    out=A-in;
    z=size(S.row,1)-in;

    %  [node,elem,face1]=v2m(b,0,10,100);
    %  exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\iceball.obj', [node(:,1)+min(S.X1) node(:,2)+min(S.X2) node(:,3)+min(S.X3)], face1)
    %  plotmesh([node(:,1)+min(S.X1) node(:,2)+min(S.X2) node(:,3)+min(S.X3)],face1)
    q=find(BW==1);
    [row col I3] = ind2sub(size(b), q);
    title(fbest, 'fontsize', 10);
    scatter3(row+min(S.x)-S.boundx,col+min(S.y)-S.boundy, I3+min(S.z)-S.boundz,'filled')
    if(~isempty(S.vessels))
    scatter3(S.row_vessel-1+min(S.x)-S.boundx,S.col_vessel-1+min(S.y)-S.boundy, S.I3_vessel-1+min(S.z)-S.boundz,'filled')
    exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\vessel.obj', [S.row_vessel-1+min(S.x)-S.boundx S.col_vessel-1+min(S.y)-S.boundy S.I3_vessel-1+min(S.z)-S.boundz], []);
    end
    hold off
    exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\iceball.obj', [row+min(S.x)-S.boundx col+min(S.y)-S.boundy I3+min(S.z)-S.boundz], []);
end
load Accumulated_Results
Accumulated_Results {end+1,1} = S.method;
Accumulated_Results {end,2} = S.bioheat;
Accumulated_Results {end,3} = S.time;
Accumulated_Results {end,4} = S.temp;
Accumulated_Results {end,5} = S.needle_num;
Accumulated_Results {end,6} = S.xr;
Accumulated_Results {end,7} = S.count;
Accumulated_Results {end,8} = [fbest z];
Accumulated_Results {end,9} = S.Total_time;
Accumulated_Results {end,10} = v;
Accumulated_Results {end,11} = S.tumor_type;
Accumulated_Results {end,12} = S.vessels;
Accumulated_Results {end,13} = S.b;
Accumulated_Results {end,14} = S.w;
save('Accumulated_Results','Accumulated_Results')